//==============================================================
// Name        : ReSEQ_bits.cpp
// Author      : Yangho Chen at University of southern California.
// Version     :
// Copyright   : Open source
// Description : C++, Ansi-style
//============================================================================
#include "ReadsMapping.h"
#include "Genome_Index_TableQ.h"
#include "ReadInBitsSet.h"
#include "ReadInBits.h"
#include "ParameterList.h"
#include "Filename.h"
#include "chdir.h"
#include "stdafx.h"
#include <string>
#include <vector>
#include <iostream>
#ifdef __GNUC__
#include <omp.h>
#endif 
using namespace std;

time_t startt, endt;
int gotoWorkingDirectory(const char* WORK_DIR) {
	if(WORK_DIR[0] == '\0') {
	    cout << "Stay in the current directory" << endl;
	} else {
		tosomedir(WORK_DIR);
	}
	Print_Working_Directory();
	return(0);
}


ParameterList getParameterList(int argc, const char** argv)
{
    // If there are none parameter   (argc == 1), load the default config file.
    // If there are one  parameter   (argc == 2), it should be the config file.
    // If there are two parameters   (argc == 3), the 1st should be RefFile, the 2nd should be Reads file List
    // If there are three parameters (argc == 4), the 3rd one should be the substitution threshold

    ParameterList parameters;
    const char* DEFAULT_CONFIG_FILE = "config.pref";

    if (argc < 2)
	{
        // Load default Setting from default file
        cout << "Try to load default parameter setting !" << endl;
        if(! parameters.load(DEFAULT_CONFIG_FILE, true)) {
        // If defatult setting doesn't exisit, ask user to input.
		    parameters.prompt();
            // Savie to the config file so user don't need to input again next time.
		    parameters.save(DEFAULT_CONFIG_FILE);
        }
	}  else if(argc == 2) {
        parameters.load(argv[1], true);
    }  else if (argc >= 3) {
        strcpy(parameters.genomeFile,argv[1]);
        strcpy(parameters.fragmentFile,argv[2]);
        if(argc >= 4) { // if argv not valid, exact match only.
            parameters.mismatchTolerance = atoi(argv[3]);
        }
        if(argc >= 5) { // if argv not valid, collect mapping statistics only.
            parameters.outputOption = atoi(argv[4]);
        }
        if(argc >= 6) { // if argv not valid, the program won't do the mapping.
            parameters.iterativeTimes = atoi(argv[5]);
        }
        if(argc >= 7) {
            strcpy(parameters.qualityFile, argv[6]); //If no valid qscore file, use no quality score for mapping
        }
    }
    return(parameters);
}

void generateOutputFileName(ParameterList& parameters ,int iterative, int run, char* filename, bool bUseQualityScore)
{
   char useQScore = bUseQualityScore ? 'Y': 'N';
   if(parameters.outputOption != 4 /* PRINT_ALL_ALIGNMENTS_OPT */)
   {
        strcpy(filename, parameters.genomeFile);
        chExtName(filename, ""); //remove extendfile
        sprintf(filename, "%s_%d_%c_%d_%d_%d.snp",\
            filename, parameters.outputOption, useQScore, parameters.mismatchTolerance, iterative, run);
   } else {
       strcpy(filename, parameters.outputFile);
   }
}

vector<string> getReadSetsFilenames(char* readSetListFilename)
{
    vector<string> readSetList;

    if(fileExist(readSetListFilename)) {
        if(hasTheExtName(readSetListFilename, ".fasta") ||
           hasTheExtName(readSetListFilename, ".fa") ||
           hasTheExtName(readSetListFilename, ".csfasta")) {
               string theOnlyReadSet( readSetListFilename );
               readSetList.push_back(theOnlyReadSet);
           return(readSetList);
        } else /* if (hasTheExtName(readSetListFilename, ".txt")) */ {
            char readsFile[FILENAME_MAX];
            ifstream readsFileList(readSetListFilename);
            while(GetNextFilenameFromListFile(readsFileList, readsFile))
            {         
                readSetList.push_back(readsFile);
            }
        }
    } else {
        cout << "Can't open the read set files list" << endl;
    }
    return(readSetList);
}

// Given a read set list and the index table, this function maps reads parallelly 
int parallelMapping(vector<string>& readSetsList,  CGenome_Index_TableQ& indexTable) 
{
    int numberOfCPUs = 1; 
#ifdef _OPENMP // Parallelization 
    numberOfCPUs = omp_get_num_procs();
#endif
    std::cout << "There are " << numberOfCPUs << " CPUs\n";
    
#ifdef _OPENMP // Parallelization 
     #pragma omp parallel for 
#endif
    for(int i = 0; i < (int)readSetsList.size(); i++) {
        CReadsMapping mapping;
        CReadInBitsSet readSet(readSetsList.at(i).c_str(), indexTable.uiRead_Length);
        TIME_INFO(mapping.mapReads(readSet, indexTable), "Mapping ");
    }
    return(numberOfCPUs);
}

/*
 * The main function construct or read in the index of a set of reference genome.
 * It takes parameters as getParameterLists();
 * If the parameter is missing, default config.pref will be search.
 */
int main(int argc, const char* argv[])
{
    const bool isMakeTable = false;	
    const bool bSaveTable = true;


	cout << "PerM program - Preprocessed Genome Beta Version!" << endl;
    gotoWorkingDirectory(""); // use default working directory

    // (1) Get parameters
    ParameterList parameters = getParameterList(argc, argv);    
    vector<string> readSetsList = getReadSetsFilenames(parameters.fragmentFile);    
    const char* firstReadSetFileName = readSetsList.at(0).c_str();
    bool bMappedSOLiDRead = hasTheExtName(firstReadSetFileName, ".csfasta");
    unsigned int readLength = get_ReadLength(firstReadSetFileName);    
    cout << "Mapping " << (int)readSetsList.size() << " read sets with read length " << (int)readLength << endl;

    // (2) Build index table
    CGenome_Index_TableQ indexTable(parameters.genomeFile);           
    if(isMakeTable) {
        indexTable.construct_IndexTable(readLength, parameters.mismatchTolerance, bMappedSOLiDRead, bSaveTable);
    } else {
    	bool sucessfullyReadTable =
            indexTable.readIn_GenomeIndexTable(readLength, parameters.mismatchTolerance, bMappedSOLiDRead);
        if(!sucessfullyReadTable) {
          indexTable.construct_IndexTable(readLength, parameters.mismatchTolerance, bMappedSOLiDRead, bSaveTable);
        }
    }

    // (3) Mapped reads
    // testGenome_Index_TableQ(&indexTable);        
    TIME_INFO(parallelMapping(readSetsList, indexTable), "");
    return 0;
}

